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Abstract 

The theoretical justification of the Hybrid Monte Carlo algorithm depends upon the 
molecular dynamics trajectories within it being exactly reversible. If computations were 
carried out with exact arithmetic then it would be easy to ensure such reversibility, but 
the use of approximate floating point arithmetic inevitably introduces violations of re- 
versibility. In the absence of evidence to the contrary, we are usually prepared to accept 
that such rounding errors can be made small enough to be innocuous, but in certain 
circumstances they are exponentially amplified and lead to blatantly erroneous results. 
We show that there are two types of instability of the molecular dynamics trajecto- 
ries which lead to this behavior, instabilities due to insufficiently accurate numerical 
integration of Hamilton's equations, and intrinsic chaos in the underlying continuous 
fictitious time equations of motion themselves. We analyze the former for free field 
theory, and show that it is essentially a finite volume effect. For the latter we propose 
a hypothesis as to how the Liapunov exponent describing the chaotic behavior of the 
fictitious time equations of motion for an asymptotically free quantum field theory be- 
haves as the system is taken to its continuum limit, and explain why this means that 
instabilities in molecular dynamics trajectories are not a significant problem for Hybrid 
Monte Carlo computations. We present data for pure SU (3) gauge theory and for QCD 
with dynamical fermions on small lattices to illustrate and confirm some of our results. 



1 Introduction 

The goal of this paper is to study the effect of errors introduced by working with finite- 
precision arithmetic on Hybrid Monte Carlo (HMC) and related algorithms. 

The HMC algorithm [l], |^, |3| allows us to generate an ensemble of field configurations 
which are selected from some probability distribution. The algorithm is "exact," in the sense 
that it is a Markov process which converges to the desired distribution with no systematic 
errors provided that the fictitious time molecular dynamics (MD) trajectories within it are 
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exactly reversible and area preserving, that the computation of the action for the Metropolis 
accept/reject step is exact, and that we have a supply of perfectly random numbers. It 
does not require that the MD integration or any conjugate gradient (CG) "inversion" of the 
fermionic kernel be carried out exactlyQ Indeed, if we were to carry out all computations 
using exact arithmetic (i.e., with no rounding errors whatsoever) then the leapfrog integra- 
tion scheme provides a method of integrating Hamilton's equations while maintaining exact 
reversibility and area preservation. 

All numerical computations carried out using floating point arithmetic are subject to 
rounding errors, but unless these errors are amplified exponentially we do not normally 
consider them to be a serious problem. The reason for this is that if a sensible rounding 
mode| is used for the elementary floating point operations then the error after a sequence 
of ./V operations, each of which gives rise to an error of magnitude e, is ey/N. The cost of 
working with d fractional digits grows approximately linearly^] with d and gives e = 10~ d . 
Therefore if we require an answer with an error smaller in magnitude than 5 in magnitude 
the cost is proportional to — \og 10 (5 / y/N) , which grows only like IniV. This is a very small 
correction^ to the growth of the cost of the HMC algorithm as the volume and correlation 
length of the system is increased. 

All the arguments we present as to the conditions under which reversibility is maintained 
apply to area preservation as well, but we shall concentrate on the former as the latter is 
less easy to determine empirically. 

2 A Brief Survey of Algorithms 

Let us very briefly summarize the various advantages and disadvantages of several variants of 
the HMC algorithm for dynamical fermion Quantum Chromodynamics (QCD) computations. 

2.1 Hybrid Monte Carlo (HMC) 

The algorithm is exact even if only a cheap approximate CG solution is carried out for com- 
puting the force along a MD trajectory, and such a CG solution is needed once per integration 
step. Nevertheless, an accurate CG solution is needed for the Metropolis accept/reject test 
once per trajectory. 

If we believe that the behavior of free field theory carries over to interacting theories like 
QCD, then choosing trajectory lengths which grow linearly with the correlation length of 

1 The initial CG vector must be chosen in a time-symmetric way. 

2 If we were to truncate all floating point numbers towards zero after each operation then the error after 
N operations would be eJV, and the cost would grow twice as rapidly with N but still logarithmically. 

3 On a typical computer chip doubling the number of digits requires doubling the number of gates in the 
adder and multiplier and increasing the clock period to allow for the worst-case carry propagation. Thus 
asymptotically the number of gates, the clock period, and the amount of memory would increase linearly 
with the number of digits d. In practice it is better to reduce the carry propagation time to grow as In d at 
the cost of increasing the number of gates as dlnd. This leads to an asymptotic cost which grows as d(\nd) 2 . 

For real computers the number of digits available is often "quantized" in units of words, and increasing 
the precision from about 7 digits (IEEE single precision) to about 14 digits (IEEE double precision) may be 
large, so asymptotic estimates must be treated with caution. 
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the system r oc £ minimizes the cost per independent configuration and give a dynamical 
critical exponent z « 1. The coefficient of proportionality depends upon which operator is 
being optimized, but it will typically be of order unity 

2.2 Second Order Langevin Monte Carlo (L2MC) 

This algorithm was introduced by Horowitz M (and is sometimes also called Kramers algo- 
rithm). It is similar to HMC except that only single step (Langevin) trajectories are taken, 
and that instead of choosing the fictitious momenta from a Gaussian heatbath after each step 
the final momenta of a trajectory are reversed and mixed with a small amount of Gaussian 
noise to serve as the initial momenta for the next trajectory. 

The advantages are that because the trajectories are short they are not subject to any 
large exponential amplification of rounding errors (which will be discussed in this paper), and 
that the the free field theory mean acceptance rate per trajectory is erfc (cVVSt 6 ) whereas 
for HMC it is eric(d VVSr^) , where V is the number of lattice sites, St is the leapfrog 
integration step size, and c and d are constants of order unity [|5], §, §• 

The disadvantages are that an accurate CG solution is required for every integration step 
(because each such step has an accept/reject test following it), and that upon a rejection the 
trajectory reverses itself (because the momenta are flipped). This last observation means 
that in order to minimize autocorrelations the acceptance rate has to be very close to one 
with a concomitantly tiny St. 

A free field theory analysis |7], §] indicates that L2MC and HMC have essentially the same 
scaling behavior in both the thermodynamic and continuum limits, so the choice between 
them must be made on the grounds of the cost of specific implementations and on the possible 
costs associated with instabilities occurring for long trajectories in HMC. 

2.3 Generalized Hybrid Monte Carlo (GHMC) 

The Generalized Hybrid Monte Carlo (GHMC) algorithm ||, || combines the HMC and 
L2MC methods giving an algorithm with three tunable parameters: the trajectory length 
r, the mixing angle between the old fictitious momenta and Gaussian noise 9, and the 
integration step size St. A free field theory analysis indicates that if St is chosen to given 
a good acceptance rate then there is a valley in the (r, 9) plane for which the dynamical 
critical exponent z = 1 and furthermore for which the cost per independent configuration 
does not vary much. 

As HMC and L2MC are both special cases of GHMC, the optimal choice of parameters 
may well not correspond to either case but in fact be somewhere in the middle of this valley. 



3 Causes of Irreversibility 



Empirical studies |H], |TT|] have shown that HMC computations for gauge theories are not 
exactly reversible. Given that they are performed using floating point arithmetic this is not 
at all surprising. However, if St or, for dynamical fermion calculations, the CG residual 



3 



are taken to be too large then very large violations of reversibility are observed. This large 
violation of reversibility was observed even when the initial CG vector was chosen to be zero. 

We recall that the HMC algorithm is exactly reversible unless the initial CG vector is 
chosen in a time asymmetric wayQ except for the effects of finite-precision arithmetic. We 
are thus presented with the question of why the rounding errors are amplified by a large 
factor. 

Three possible causes of this amplification suggest themselves: 

• A small error changes the Metropolis accept/reject choice; 

• A small error changes the number of CG iterations; 

• The MD trajectories exhibit an exponential sensitivity to initial conditions (in other 
words they are unstable or chaotic) [|ll|, 0. 

It is easy to see that although the first two mechanisms produce a large change they oc- 
cur correspondingly infrequently, and therefore are not compatible with the nature of the 
irreversibility observed. We therefore conclude that the irreversibility observed is caused by 
instabilities in the MD trajectories. 

One may wonder to what extend irreversibility introduces systematic errors into the 
results of an HMC calculation. If the irreversibility is due to some underlying chaos in the 
equations of motion perhaps this just has the same effect as a different choice of random 
numbers for the fictitious momentum heatbath? Indeed, it is hard to see significant effects 
on physical observables unless truly huge violations of reversibility are induced (in which case 
the trajectories often become so unstable as to cause numerical overflows). We did observe^] 
that when the trajectories become irreversible there is often a large change in the fictitious 
kinetic energy as well (even though the total fictitious energy SH stays small). This means 
that the distribution of fictitious momenta is clearly wrong, and thus the computation is 
erroneous even if the effect on physically interesting observables is hard to detect. 

4 Causes of Instability 

There are two distinct mechanisms which cause the MD equations to become unstable: 

• The discrete integration procedure (leapfrog) diverges exponentially from solutions of 
the classical equations of motion. This instability should grow with the number of 
integration steps, and thus have a short characteristic time scale for large volumes. 

5 If the CG solution was found exactly then it would of course be independent of the initial vector. In 
practice only enough CG iterations are carried out to reduce the residual to below some preset value, and 
the Krylov space thus explored depends on the initial vector. If this initial vector depends on the solution 
for the previous step then these Krylov spaces will be different for forward and backward trajectories, and 
we will induce reversibility errors depending on the CG residual used. 

6 See for example Figure || on page [l5| 
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• The classical equations of motion are themselves chaoticQ. These may be expected to 
have a time scale independent of the step size, and presumably determined by some 
characteristic length scale of the system. 

In section |5] we analyze the instabilities of the leapfrog integration scheme for free field 
theory, and in section ^] we consider the chaotic nature of the underlying continuous time 
equations of motion themselves and propose a hypothesis as to how the time scale of the 
corresponding instabilities behaves as the system approaches the continuum limit. 



5 Free Field Theory Analysis 

The phenomenon of the leapfrog integration scheme diverging from the true trajectory occurs 
even for free field theory, and it is instructive to see what happens there f7|, || . 

Consider a system of harmonic oscillators {<p p } for p e Zy. The Hamiltonian on fictitious 
phase space is 

5.1 Single Mode Analysis 

The Hamiltonian is diagonal, so we may temporarily consider the evolution of a single mode 
and set its frequency u p = 1. In fact, we will find that the essential features of the instability 
are displayed by a single mode. 

The leapfrog discretization of Hamilton's equations can be written as the following matrix 
acting on the phase space vector ((f), n) 



U{5r) 



1 - \5t 2 5t 
-St + ±5r 3 1 - \5t 2 



The general area-preserving reversible linear map on ((f), 7r) phase spacef] may be parameter- 
ized as 

U(5t) = ( C0S K^)5t] am[ p$) 5T] 
\ — p(5r) sm[K,(5r)5r\ cos[k(8t)8t] 

where k and p are even functions, but not necessarily real valued. For the lowest-order 
leapfrog algorithm k and p are 



. . . cos (1 — \8. , . . . , 



7 Chaos in classical Yang-Mills dynamics (but not in fictitious time) has been investigated by several 
authors. See and references cited therein. 

8 If F = ^ J is the momentum-flip operation then F~ 1 U(St)F — U(—St), from which it im- 

mediately follows that Ui i and Z7 2i 2 are even functions of St, and U\ t 2 and 1/2,1 are odd functions of St. 
Furthermore U(St)U{—St) = I implies that U\ t \ = 1/2,2 except for the trivial case where C/1,2 = C/2,1 = 0. 
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With this parameterization it is easy to see that the evolution operator for a trajectory 
of t/St leapfrog steps is 



U(r) 



r /r \ i sm[K( 5t)t] 
cos[k{St)t] 

-p(oY) sin(/t(<5T)r] cos[k(St)t] 



and thus the change in fictitious energy over this trajectory is[] 

/ , \ T 



SH 



1 



E 

peZ v 



7T-, 



u?u p 



I 



(pp 



Tin 



It is clear that some kind of critical behavior will occur at St = 2. 

Given this evolution matrix we can compute the probability distribution of SH, 



Psh(0 



1 

z 



[d(f)][diT]e~ H 5(t-5H) 



drj 



[d(f)][d7r]e- H 
dr) 



dr) 
-oo 2tc 



where T 



— 27r U P y/detp + iri(UjU p - I)] J-°° 2vr Up ^ + 7^(1 - ^) 
(p — ~) sin(«r) • It is easy to see that T rapidly approaches its asymptotic 



behavior of growing proportionally to e UT where 



±— Rein 

or 



\St 2 - 1 ± SrJiSr 



(2) 



(which falls as 21n(oY)/oY as St — > oo). For <5r < 2 not only is the exponent v zero, but T 
is bounded by a constant. That the exponent decreases as St — > oo merely reflects the fact 
that T is growing exponentially in the number of MD steps but only algebraically in the step 
size St. This is illustrated in Figure [TJ 

For a single mode the probability distribution of SH may be evaluated by changing 
variable to r)' = r\ + i/2 



Psh(0 



drj 



271 Jl + T 2 irj(l - irf) 



Mi 



°° +i /2 drj' 

i/2 IT 



cos(r/£) 



(3) 



The singularities in the integrand occur when 1 + T 2 {\ + r/' 2 ) = 0, for which r/' 2 < — \ and 
hence | Im?/| > \. We may thus shift the integration contour to the real axis, and defining 



x = Trj'/Jl + \T 2 we obtain 



dx 



7tT Jo VT+~? 



cos 



xt; 

y 



4 

J>2 



which can be expressed in closed form in terms of a modified Bessel function 



Psh(0 



ttT 
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Figure 1: The "exponent" Reln(T)/r is shown as a function of the integration step size St 
and the trajectory length r. This quantity approaches v as r — > oo. 




dH 



Figure 2: The logarithm of the probability distribution for the change in fictitious energy 
SH for a single mode is shown as a function of SH and the MD leapfrog integration step size 
St for a trajectory of length r = 4 (rounded to the nearest integer multiple of St). 
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Figure 3: Acceptance rate for a single mode as a function of the MD step size St and 
trajectory length r (rounded to the nearest integer multiple of St). 



The logarithm of the acceptance rate is shown in Figure ^|. 
The average Metropolis acceptance rate is 



. POO 

mm(l,e" SH ))= / d£ min(l, e^)P 5H (0- (4) 



For a single mode we may use eq. ([3]), and evaluating the integral over £ gives 

2 r 00 dx 



p _ 

ttT 7-oo v/fT^ [l + x 2 (l + ^r) 



setting x = | — an d f = ^ an f we easn y obtain 



sin# r°° dy 6 2 , / 2 

' ■ tan ' 



7T Jo y 2 + 2ycos^+l 71 71 \T 

and this is shown in Figure ^ as a function of the step size and trajectory length. Observe that 
there is a "wall" at St = 2 which separates the region where the acceptance rate oscillates 
as a function of r from the region where it plummets exponentially. 



5.2 Multiple Mode Analysis 

These results can be generalized to a field theory with many "stable" modes and some number 
of "unstable" ones. In the general case let us consider the acceptance rate P acc defined in 
eq. and use the expression of eq. (P for the distribution of values of SH. As before we 

9 The subscript p indicates that all fictitious times are to be multiplied by the frequency u> p . 
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may make the change of variable 7/ = 77 + i/2 and observe that there are no singularities in 
the strip of the complex 7/ plane satisfying | Imr/| < \. We thus obtain 



/0 roo 
-oo Jo 

drf 1 



- 2vr n p ^i + Tp 2 (i + r/2) 

30 drj' 1 



00 



27r (i + v 2 )n P vi + T p 2 (i + ^ 



Let us now assume that there are a small number N of "unstable" modes with T q 3> 1 
and a large number V of "stable" modes with T p C 1 such thatQ Y%=i Tp = er~ 2 and 
Z)^=i Tp" — 0(V 1 ~ n ). We can thus compute an asymptotic expansion valid for large V, 



P„, 



dj ^p {-§ E p y = i In [l + T p 2 (i + r/' 2 )] } 
00 2vr ( 1 + ^2) jjtf /l + T^i + r/' 2 ) 



e -V** 3 / ~ ^_ t fuiV tVi _i_ V 2 ^ 2 

00 2vr ( 1 + n<? ^l + T^i + r/' 2 ) 



i + iE^ 4 a+v 2 ) 2 + 



where the "stable" modes correspond to subscript p and the unstable ones to subscript q. 
If there are a few modes with V~ l / A <C uj p 5t <C 1 then these may be taken into account by 
computing further terms in the asymptotic expansion above. On the other hand, we may 
carry out an expansion in powers of 1/T q for the "unstable" modes by writing 

1 (1 + „/2-S -iV/2 / 1 \ 

n-,/ 1 + ^ + ,< k l^( 1 -^^ + " rl + -> (5) 



Using this expansion in the integrand above reduces both the large T and small T expansions 
to a sum of integrals of the form 

jTdi/e-^/^Ci + O-* 72 = 2 fc ~ 2 v^, ^) (6) 

where [/ is a confluent hypergeometric function.^ It is easily verified that for N = the 
leading term in the asymptotic expansion is just erfc(l/v / 8a 2 ). 

In Figure |] we show the P acc for the case where there are many "stable" modes with 
a = 1 as a function of the value of a for one additional "unstable" mode. It is clear that 
the expansions calculated above are very good as long as T is not close to one. In Figure |5] 
we plot In P acc as a function of T for various values of N, the number of "unstable" modes. 
Each additional unstable mode clearly reduces the acceptance rate by a large factor for these 
parameters, and this makes it reasonable to assume that at least the qualitative behavior of 
the instabilities in the system is apparent from the case where there is just a single unstable 
mode. In Figure ^| we show a comparison of the dependence of the mean acceptance rate on 

10 If St < 1 then T p = -\{lo p 8t) 2 sin^r) + 0(5t 4 ), and thus ^ =1 T p 2 = kVSr 4 + 0(8t 6 ) where the 
spectral average k = j^y JZp—i oj p sin 2 (uj p t) is of order unity. In order to obtain a non- negligible acceptance 
rate we must take St to be of order V~i, in which case the quantity a will also be of order one. 

"This is related to Whittaker's function W\^(z) = z» + ie~ z / 2 U(n - A + ±, 2/i + 1; z). 
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Figure 4: Acceptance rate for many "stable" modes with a = 1 and one "unstable" mode as 
a function of the value of T for the "unstable" mode. In addition to the result obtained by 
numerical integration various orders of the small- and large-T expansions are shown. 



Log Pace 




Figure 5: Logarithm of the mean acceptance rate for N "unstable" modes as a function of 
a for T = 25.63 which corresponds to 5t = 2.1 and r = 9. The four curves are for N from 
(top) to 4 (bottom). 
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Figure 6: Comparison of the acceptance rate for a single mode (multiplied by a factor of 
erfc(l/v / 8<j 2 ) compared with the corresponding acceptance rate for one "unstable" mode 
and many "stable" modes for o = 0.5, 1, and 1.5. 

T for a single mode and for a system of many "stable" modes and one "unstable" one. This 
demonstrates that the behavior of a single mode dominates the instabilities in the case of 
free field theory. 

We notice that in order to keep the acceptance rate constant as the lattice volume V — » oo 
we must decrease 6r so as to keep V5t 4 fixed. Thus as we approach the thermodynamic 
limit the instabilities will go away: in other words the leapfrog instability is a finite volume 
effect. 

5.3 Interacting Field Theories 

When interactions are present it is no longer meaningful to talk of independent modes, but 
for weak interactions like those for asymptotically free field theories at short distance it is 
still useful to consider the system as a set of weakly coupled modes. We expect that the 
onset of the integration instability will still be caused by the highest frequency mode, and 
will occur at u max 5r = 2. The forces acting on the highest frequency mode due to the other 
modes will fluctuate in some complicated way, and consequently we would expect that the 
"wall" at uj max 5r = 2 gets smeared out. 

In order to illustrate this effect we have made numerical studies of a simple model of a 
single harmonic oscillator whose frequency is randomly chosen from a Gaussian distribution 
with mean unity and standard deviation o before each MD step. The behavior of this model 
is just like that of the free field theory considered above, except that the "wall" at 5r = 2 
gets spread out. Our numerical results are shown in Figure [?|, where the analytic result for 
a = is taken from eq. (|2]). We conjecture that interacting asymptotically free field theories 
like non-Abelian gauge theories behave similarly, and this is supported by our numerical 
results. 
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Figure 7: Results from fluctuating-frequency harmonic oscillator model. The curve for a = 
is given by eq. (H) on page §. 
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We do not expect the introduction of dynamical fermions to produce any qualitatively 
different behavior. For QCD computations using the HMC algorithm the pseudofermions 
produce a force which we expect to be proportional to l/mj, so the "highest frequency" of 
the system, which is responsible for the integration instability, will grow as oo max oc 1/rrif 
too. In particular we expect u max to become large as rrif — > 0, so the critical value of St will 
become small for Wilson fermions. An inaccurate CG solution will produce errors 

which are multiplied by the scale of this fermionic force, and will become important when 
the size of the errors in the CG solution become comparable with the size of the fluctuations 
in the forces due to the fields themselves. 

It is interesting to observe that reversibility can be violated significantly even though 
the acceptance rate does not become small. For example we will see an example of this in 
Figure [TO on page ITT] for St = 0.005 and r = 50. This is because the leapfrog integration 



scheme tries to conserve energy even when it has wandered far away from the true continuous 
time trajectory. It is therefore prudent to verify that reversibility is satisfied to the precision 
required, even when the HMC algorithm has an average acceptance rate which is close to 
one. 



6 Chaos in Continuous Time Evolution 



Jansen and Liu |TTJ observed that there is a second cause of instability in the MD trajectories, 
which occurs even when the integration step size St — > 0. This is because the underlying 
continuous fictitious time equations of motion are themselves chaotic. In the chaotic regime 
two nearby classical trajectories diverge exponentially with a characteristic exponent v called 
the Liapunov exponent. Obviously such behavior cannot be studied in the context of free 
field theory. 

The instability due to chaos cannot be removed by reducing St, and is therefore not a 
finite size effect. If the system we are studying exhibits critical slowing down,0 and we 
believe that the mechanism for reducing this to z = 1 in free field theory is applicable to this 
system too, and we wish to use the HMC algorithm, then we must scale the trajectory length 
with the correlation length, r oc £. This means that there is the potential for exponential 
amplification of rounding errors. 

Even if there is such an exponential amplification of rounding errors, reversibility can 
be maintained to any desired precision by a linear increase in the number of digits used 
for floating point arithmetic. This is no longer a negligible addition to the computational 
complexity of the problem (i.e., it is not just a logarithmic correction), but it is still a small 
correction to the power dependence of the computational cost on the correlation length and 
volume. 

This problem can be avoided using the GHMC algorithm, but not only is this unnecessary 
until chaotic amplification of rounding errors becomes important, but also as we shall see it 
might never be necessary at all. 

The reason for this is that the Liapunov exponent (averaged over the equilibrium distri- 
bution) is not constant as a function of (3. If the chaotic dynamics is not only a property of 
the underlying continuous fictitious time evolution, but is also a property of the underlying 



2 This is probably not the case for most current lattice QCD computations. 
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(space-time) continuum field theory, then the Liapunov exponent would be constant when 
measured in "physical" units, that is z/£ would be constant as £ — > oo. In other words this 
hypothesis would say that v oc l/£, or for QCD with rif flavors of fermions 

uoce~^ with A) = ^4^, (7) 

where the last relation follows from perturbation theory for f3 large enough that the system 
exhibits asymptotic scaling. In section [7.2| we shall present numerical evidence that this 
hypothesis may indeed be true. If this is the case, then tuning the HMC algorithm by 
varying the trajectory length proportionally to the correlation length does not lead to any 
change in the amplification of rounding errors as we change £. 



7 Monte Carlo Results 

We have carried out extensive numerical studies of reversibility errors for (quenched) SU(3) 
gauge theory and for full QCD with two flavors of dynamical Wilson fermions. Measurements 
were made principally on 4 4 lattices, but additional data on 8 4 lattices was used for an 
analysis of finite size effects. 



7.1 Numerical Determination of Reversibility Errors 

In order to measure the deviation from reversibility of MD trajectories we chose some initial 
point in fictitious phase space (Ui, m) and followed an MD trajectory of length r to the final 
point (Uf, TTf); we then reversed the fictitious momenta and followed the backward trajectory 
of length t from (Uf, —TTf) to (U-, 7r-). At the end of this backward trajectory we measured 
the amount by which it failed to return to its starting point by computing two quantities, 
the change in energy ASH = H(U[, 7r|) — H(Ui, 7Tj), and the norm of the change in the gauge 

fieldF] \\ASU\\ = \\U! - Ut\\, where \\U\\ 2 = Ex^a^P^Uf , a and b being SU(3) color 
indices. Observe that both quantities are extensive in the lattice volume, that A5H is gauge 
invariant whereas ||A5C/|| is not, and that ||A5C/|| cannot be small unless ASH is too. These 
quantities can also be considered as the difference between the change in energy or gauge 
field over the forward and the backward trajectories, which serves to explain our notation. 
Both of these quantities vanish identically for reversible trajectories, of course. 



7.1.1 SU(3) Gauge Theory 

For SU(3) gauge theory an example of a long trajectory is shown in Figure ^. This shows 
how the energy SH changes along a long trajectory. The change in energy over the forward 
trajectory was SH = 1.2 (corresponding to an acceptance probability of 30%), which differed 
from the change in energy over the backward trajectory by ASH = 0.36. This clearly 

13 The configuration norm ||A5J7|| could have been made into a phase space norm by adding ||A#7r|| to it, 
but we do not believe that this would make any significant difference in practice (and it would require saving 
the initial momenta iri , which are not needed in the usual HMC algorithm but are needed for the L2MC and 
GHMC algorithms). 
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Figure 8: The contributions to SH from the fictitious kinetic energy and the potential 
energy (gauge action) for a trajectory of length r = 40 with step size St — 0.1 at (3 — 
5.7. The two graphs show the contributions to the energy for the first and last four units 
of fictitious time respectively, overlaid with the corresponding quantities on the backward 
trajectory. The graph on the right shows that the forward are backward trajectories are 
essentially indistinguishable immediately after the reversal, and the graph on the left shows 
how the backward trajectory finally deviates from the forward one. Note that the cancellation 
between SKE and SPE is still at the 2% level even at the end of this very long trajectory 
with large step size, as the SH values have been multiplied by a factor of 100. 

shows that the leapfrog integration scheme tries very hard to conserve energy even though 
it has wandered far from the true path: indeed, the backward trajectory ends up a distance 
||A<5£/|| = 20 away from the starting configuration. This example demonstrates that there 
can be severe violations of reversibility even though the acceptance rate is quite reasonable, 
but it should be noted that a trajectory length of r = 40 is unreasonably long for a system 
whose correlation length £ = 0(1), especially on a 4 4 lattice! 

By carrying out such forward and backward trajectories for a variety of step sizes and 
trajectory lengths we produced Figures || and [10]. In the former we plot the logarithm of 
\A5H\, \5H\, and ||A5C/|| as a function of the integration step size. The differences between 
measurements made on different equilibrated configurations were very small. The "kinks" 
in some of the curves are just a consequence of rounding r to the nearest integer multiple of 
St. The top and bottom graphs clearly show the integration instability "wall" at St ~ 0.6, 
which has spread out just as we would expect. The middle graph shows that by the time 
one has reached the "wall" SH = O(10 3 ), so the integration instabilities are of no practical 
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Figure 9: Results for pure SU(3) gauge theory with (3 = 5.7 on a 4 4 lattice as a function of 
5r. 



importance for this system. If one looks carefully at the left end of the bottom graph one 
sees that ||A5C/|| increases a little as St — > 0. This is because the number of integration 
steps n = t/St is increasing, and therefore rounding errors are being amplified by the usual 
sjn factor. 

For extremely long trajectories reversibility is violated even for very small values of St; 
in order to understand this it is instructive to look at Figure [TO, where we plot the logarithm 
of SH and ||A5Z7|| as a function of the trajectory length for two values of (3 and three values 
of St. The lower graph immediately shows us that even for very small step sizes the errors 
are amplified exponentially, with an exponent which is independent of St (i.e., the lines for 
the same (3 value but different St are parallel). This is evidence for chaos in the continuous 
time equations of motion, as suggested by Jansen and Liu |TTJ. The lines for small step 
size on the lower graph eventually curve downwards for the largest r values because SU (3) 
is a compact manifold so there is a maximum distance two configurations can be from one 
another, and this bound is becoming saturated. The upper graph shows that SH stays small 
for small St even though the rounding errors have been enormously amplified. 



7.1.2 QCD with Wilson Fermions 

We have carried out the corresponding analysis for full QCD with two flavors of Wilson 
fermions with configurations taken from two ensembles, one at (3 = 5.1 and k = 0.16 which 
corresponds to quite heavy quarks, and one at (3 — 3.5 and k = 0.2225 which is extremely 
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Figure 10: Results for pure 577(3) gauge theory with j3 
function of r. 



5.7 and 6.2 on a 4 4 lattice as 



close to k c . We took care that both ensembles were in the "confined phase," if we use ther- 
modynamic terminology to approximately describe the finite size effects on our symmetric 
lattices. 

Figure [II] shows a typical trajectory for the light dynamical fermion system. For this 
trajectory 5H = 3.9 x 1CT 2 corresponding to an acceptance probability of 96%, and 1 1 A6U\ \ = 
1.4 x 10~ 2 . This trajectory exhibits some typical behavior: the potential energy oscillates 
smoothly with a period of about 2 MD time units; the fermionic energy is much more 
jittery; and the kinetic energy adjusts itself to cancel the other two components leaving 
a small residual 5H whose magnitude is about a percent of that of its constituent parts. 
The three constituent contributions to the energy are of similar magnitudes, which indicates 
that the fermions are indeed very light: with heavier fermions the fermionic energy is much 
smaller. 

From many such trajectories on several equilibrated configurations we produced the data 
shown in Figures |I~2| and |T3|. For simplicity^ we only show the data for the heavy fermion case 
in Figure 12. The results are very similar to the pure 577(3) gauge theory results of Figures ||] 
and [10], especially if one notes that the trajectory lengths are somewhat longer in the present 
case. Just as in that case we observe the smeared-out "wall" at which the leapfrog integration 
scheme becomes unstable, but that the acceptance rate is already completely negligible by 



14 The data for the light fermions is similar in form but on a very different St scale. The results obtained 
from all the data will be shown in Figure |l4] on page |2l]. 
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Figure 11: The contributions to 5H from the fictitious kinetic energy, potential energy (gauge 
action), and the fermionic energy for a trajectory of length r = 7 with step size 5t = 0.01 
at /3 — 3.5 and k = 0.2225. Each curve is overlaid with the corresponding quantities on 
the backward trajectory, but the forward and backward trajectories are indistinguishable on 
this scale. Note that 5H has been scaled by a factor of 100 because it has a magnitude of 
less than a percent of its constituent parts. The number of CG iterations used to reach the 
residual r = 10~ 8 from a zero start is also shown (divided by a factor of 10). 
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Figure 13: Results for QCD as a function of trajectory length r. The data shown is from 
4 4 lattices with a CG residual r = 1CT 8 . 

the time the wall is reached. We also see that long trajectories with small step sizes exhibit 
chaotic behavior. 

All our data show a clear exponential instability in 1 1 A5U\ | as a function of r, so we fitted 
the data over the range of r for which this exponential behavior was obvious and extracted 
a characteristic exponent v (which we shall henceforth call the Liapunov exponent, although 
this might not always be quite the correct terminology) . In Figure |T3| we show the Liapunov 
exponent as a function of the step size St for pure SU (3) gauge theory, and QCD with heavy 
and light dynamical fermions. The results show the same qualitative behavior as that of our 
toy model which was exhibited in Figure on page [12], except that 

• The "wall" where the integration instability sets in is at a different value of 6t (about 
0.6 for the pure gauge theory and heavy fermion cases, and about 0.1 for the light 
fermion case). This probably just reflects the different highest frequencies of the sys- 
tems. 

• The exponent is not zero for small 6t, but has some fixed non- vanishing value. This 
is just the evidence for chaotic continuous time dynamics discussed earlier. 

The integration instability is easily avoided by choosing a sufficiently small step size, which 
is forced upon us anyhow as the lattice volume becomes large if we want to maintain a 
reasonable Metropolis acceptance rate. This, of course, is exactly what we mean by saying 
that the integration instability is a finite volume effect. Furthermore, it is clear that even 
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Figure 14: The Liapunov exponent v is shown as a function of the integration step size 
5t. The top graph is for pure 577(3) gauge theory, the middle one is for QCD with heavy 
dynamical Wilson quarks, and the bottom one is for QCD with light dynamical Wilson 
quarks. Note that the scale for the light quark case is quite different from the other two. 
The error bars show the standard deviation of measurements made on three independent 
configurations. 
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Figure 15: ||A5Z7|| as a function of trajectory length r for a variety of values of (3 for pure 
SU (3) gauge theory on 8 4 lattices. 



a 4 4 lattice is sufficiently large for us to be driven well away from the "wall" for all the 
parameters we have considered. 

We shall therefore turn our attention to the chaotic instability which is present for any 
value of 5t, and in the next section we shall investigate the dependence of v upon the 
parameters of our lattice field theories. 



7.2 Parameter Dependence of Liapunov Exponent 

We investigated the dependence of the Liapunov exponent v on the coupling constant (3 for 
pure £77(3) gauge theory.^ The results are shown in Figure [15] and [TJ]. In the former we 
plot In ||A5C/|| against r for a variety of values of (3. All of the curves exhibit a very clear 
exponential instability: for small j3 the relation is almost completely linear on this semi- 
logarithmic plot, whereas for large (3 there is some initial curvature before a linear region 
is reached. The slope (i.e., the Liapunov exponent) decreases very rapidly between (3 = 5.4 
and 6.0, which is just below the finite-size analogue of the finite temperature deconfinement 
transition on a 8 4 lattice. 



In Figure [16] we plot the measured values of the Liapunov exponent v as a function of j3 



for both 4 4 and 8 4 lattices. For small (3 the lattice theory is in the strong coupling regime 

15 Somewhat similar results have been reported by Jansen and Liu Ji"^ , although their data does is not 
obviously consistent with ours, especially for large (3. 
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Figure 16: The Liapunov exponent v is shown as a function of (3 for pure SU(3) gauge theory 
on 4 4 and 8 4 lattices. The curve is a one-parameter fit of the functional form given in eq. ([?]) 
to four of the 8 4 data points (f3 = 5.4, 5.5, 5.6, and 5.7). The data was measured on three 
4 4 configurations and two 8 4 configurations. 

and does not obey the asymptotic scaling behavior mandated by the renormalization group 
equations and the perturbative /3-function. For large (3 the system is in a tiny box, and is 
thus in the deconfined phase. Therefore we can at best only trust the data in some "scaling 
window" near j3 — 5.7. In this region we have fitted the 8 4 data to our suggested asymptotic 
scaling form of eq. (0) and we find the fit surprisingly good, especially considering that there 
is only one free parameter. It is dangerous to rely too much on asymptotic scaling results 
obtained on such small lattices, so we will content ourselves with the statement that our data 
is consistent with our hypothesis that u£ is a constant. The value of v is also sufficiently 
small that the amplification of rounding errors for trajectories of length r = £ is a factor of 
0(1), and thus unimportant. 

We do not have a clear understanding of why the value of v at large (3 is so large for 4 4 
lattices, other than the obvious suggestion that this is a finite size artefact. The data from 
8 4 lattices shows a small decrease in v for large /3, but this is not convincing evidence that 
the data will eventually fall onto our hypothetical asymptotic scaling curve. 

We can gain a little more understanding of the mechanism by which v could decrease as 
we approach the continuum limit by noting that v depends on (3 in two ways: there is an 
explicit /^-dependence of the equations of motion, and there is an implicit /^-dependence in 
the equilibrium ensemble of configurations over which v is measured. This latter effect is 
illustrated in Figure [T7], where we plot In | |A5Z7| | against r for (3 = 5.1 and k = 0.16 on a 4 4 
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Figure 17: In ||A5i7|| is shown as a function of trajectory length r (with (3 = 5.1, k = 0.16, 
and 5t = 0.1 on a 4 4 lattice) for starting gauge field configurations which are "hot" (all link 
variables U^{x) chosen independently according to SU(3) Haar measure), "equilibrium" (a 
typical configuration chosen from the equilibrium ensemble), and "cold" (all link variables 
U^{x) = 1). Fits of the form ||A<5£/|| = ke UT for the hot and equilibrium configurations and 
||A5C/|| = kr 3 / 2 for the cold configuration are shown. 

lattice for three different configurations, one hot, one cold, and one chosen at random from 
the equilibrium distribution.^ The hot configuration yields a larger Liapunov exponent than 
the equilibrium one, and the cold configuration is consistent with a power-law dependence 
of ||A5Z7|| on r. This result suggests that we should not expect to get the correct continuum 
result from lattices which are too small to be in the confined phase (even if we wanted to 
consider a finite temperature system in the deconfined phase we should at least ensure that 
the spatial volume is large enough). 

In Figure [17| we also show two-parameter fits to the hot and equilibrium trajectories of the 
form || A8U\ | oc e UT , and a one-parameter fit to the cold trajectory of the form 1 1 ASU\ | oc r 3 / 2 . 
These fits are stable if we add more free parameters by fitting to a combination of an 
exponential and power dependence. The power-law behavior of the cold trajectory may 
be explained by a combination of a factor of yfr for the random walk evolution of rounding 
errors and a linear factor characteristic of the divergence of nearby trajectories for interacting 
dynamical systems in the stable (non-chaotic) region of phase space. . 

16 We verified that the results for the equilibrated configuration did not change significantly when we 
selected another such configuration. 



24 



1.4 - 




0.8 - 
-10.5 - 




-14 -12 -10 -8 -6 -4 -2 

'Oflio r 



Figure 18: The dependence of P acc and ||A<5£/|| on the CG residual r. We write ||A5Z7|| = 
ke UT , where the coefficient k and exponent u are functions of r. A time-symmetric zero initial 
CG vector was used. The data is for 4 4 lattices at (3 — 3.5 and n = 0.2225 with 5t = 0.01. 



7.3 Reversibility and Conjugate Gradient Accuracy 

Our final topic is to investigate the effect of inaccurate CG solutions on instabilities and 
reversibility for dynamical fermion computations. 

In Figure |18| we show the dependence of the Liapunov exponent v and the Metropolis 
acceptance rate as a function of the logarithm of the CG residual for a time-symmetric (zero) 
initial CG vector. We define the CG residual for the approximate solution x' of the linear 
equations Ax = b as r = \\b — We also show the logarithm of the coefficient 

k(r) for completeness. The results clearly show that the choice of residual has no effect for 
r < 10~ 5 , and that for r > 10 -2 the acceptance rate is essentially zero. There would seem to 
be no benefit from finding the solution more accurately than is needed to give a reasonable 
acceptance rate. 

The effect of improving the convergence of the CG algorithm by using a time-asymmetric 
initial vector is illustrated in Figure pl| For a residual r = 10~ 8 the solution vector is 
sufficiently independent of the starting vector that the choice of initial vector has no effect 
on the reversibility (or otherwise) of the trajectory. On the other hand, for r = 10 -5 there 
is a large difference between the time-symmetric start and the intrinsically irreversible time- 
asymmetric starts. 

It is clear that the reduction in the number of CG iterations required to reach a prescribed 
residual obtained by using an starting vector derived from previous nearby solutions |14| 
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Figure 19: The dependence of In ||A5C/|| on r for different choices of initial vector for the 
CG algorithm. The results for the time-symmetric zero start are to be compared with those 
for which the initial CG vector was chosen to be the solution from the previous time step 
(constant interpolation), and a linear extrapolation from the previous two time steps (linear 
interpolation). For a residual r = 10~ 8 the initial vector has little effect upon the solution 
obtained, but for the larger residual r = 10 -5 the explicit violation of reversibility is very 
significant. The data is from 4 4 lattices with (3 = 3.5, k = 0.2225, and 5r = 0.02. 
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needs to be offset against the need to have a more accurate solution in order to avoid 
increased reversibility errors. To reach definitive conclusions more extensive studies would 
have to be done. 



8 Conclusions 

• Instabilities and the concomitant amplification of rounding errors caused by inaccu- 
rate numerical integration of the equations of motion should not be a problem if the 
integration step size is chosen suitably. 

• The instabilities due to integration errors become unimportant as V — > oo. 

• We hypothesize that the Liapunov exponent falls as v oc £ _1 for sufficiently large cor- 
relation lengths £, and therefore we can choose trajectories of length r oc £ which may 
reduce critical slowing down without suffering from exponentially large amplification 
of rounding errors. 

• The measured values of the Liapunov exponent are of 0(1), so the exponential ampli- 
fication factor for trajectories of length r « £ are this only e"* = 0(1). 

• We have found examples with large reversibility errors but reasonable acceptance rate 
and small \ A5H\. It is always prudent to verify that ||A5Z7|| is small. 
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